library(ggplot2)These set up the phenotype and its peak location for further analysis.
chr_id <- "11"
scan_window <- c(96,98)
pheno_name <- "LD_light_pct"
peak_Mbp <- mean(scan_window)
window_Mbp <- diff(scan_window) / 2local <- TRUE # consider local QTL (from mRNA traits)
other <- FALSE # consider other types of phenotypes if TRUE
qtls <- 1 # or 2 for number of driversThe qtl2shiny package uses a data frame of project
information (folder locations).
dirpath <- "../../qtl2shinyApp"
project_info <- data.frame(taxa = "CCmouse",
project = "Recla",
directory = "qtl2shinyData")These encode a hierarchy:
datapath <- file.path(file.path(dirpath, project_info$directory))
taxapath <- file.path(datapath, project_info$taxa)
projectpath <- file.path(taxapath, project_info$project)Information for all phenotypes: all covariates and the
physical map pmap.
covar <- readRDS(file.path(projectpath, "covar.rds"))
str(covar)## 'data.frame': 261 obs. of 6 variables:
## $ sex : chr "male" "male" "male" "male" ...
## $ Cohort : chr "Fall" "Fall" "Fall" "Fall" ...
## $ Group : chr "1" "1" "1" "1" ...
## $ Subgroup : chr "1B" "1A" "1A" "1A" ...
## $ ngen : chr "4" "4" "4" "4" ...
## $ coat_color: chr "AG" "AG" "AG" "WH" ...
pmap <- readRDS(file.path(projectpath, "pmap.rds"))
sapply(pmap, range)## 1 2 3 4 5 6
## [1,] 3.252796 3.164247 3.163885 3.427907 3.160082 3.278808
## [2,] 194.886567 181.812247 159.490311 156.138564 151.428385 149.460545
## 7 8 9 10 11 12
## [1,] 3.215308 3.235539 3.582987 3.180584 3.355458 3.562561
## [2,] 145.309341 128.662187 123.970379 130.550456 121.946117 119.342990
## 13 14 15 16 17 18 19
## [1,] 4.466555 3.262695 3.269262 3.16841 3.264958 3.235759 3.286882
## [2,] 119.392941 124.681199 103.176510 97.89090 94.717592 90.441479 61.167435
## X
## [1,] 5.503248
## [2,] 169.979596
Now subset on region of interest:
pmap <- pmap[[chr_id]]These are used to find objects that are not immediately loaded, including genotype probabilities and gene and variant information. See next section on query routines.
The phenotype data are extracted here.
pheno_data <- readRDS(file.path(projectpath, "pheno_data.rds"))
dim(pheno_data)## [1] 261 26
colnames(pheno_data)## [1] "OF_distance_first4" "OF_distance"
## [3] "OF_corner_pct" "OF_periphery_pct"
## [5] "OF_immobile_pct" "OF_center_pct"
## [7] "LD_distance_light" "LD_transitions"
## [9] "LD_light_first4" "LD_light_pct"
## [11] "VC_top_distance_first4" "VC_top_time_first4"
## [13] "VC_top_distance" "VC_top_time_pct"
## [15] "VC_top_velocity" "VC_bottom_distance_first4"
## [17] "VC_bottom_time_first4" "VC_bottom_distance"
## [19] "VC_bottom_time_pct" "VC_bottom_velocity"
## [21] "VC_bottom_transitions" "HP_latency"
## [23] "TS_frequency_climbing" "TS_time_immobile"
## [25] "TS_latency_immobile" "bw"
The peaks and analyses_tbl have information
on phenotypes. The peaks has location of all previously
found peaks, used to find what maps in a region. The
analyses_tbl has details of transformations and covariates.
Phenotypes are ordered into groups and types, which is useful for
searching and for determining sets of covariates and transformations
when there are many phenotypes to consider.
peaks <- readRDS(file.path(projectpath, "peaks.rds"))
str(peaks)## 'data.frame': 450 obs. of 8 variables:
## $ pheno : chr "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" ...
## $ chr : Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ pos : num 84.1 115.9 26.5 21.7 128.5 ...
## $ lod : num 5.25 5.67 6.55 3.65 5.03 ...
## $ longname : chr "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" ...
## $ output : chr "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" "OF_distance_first4" ...
## $ pheno_group: chr "group" "group" "group" "group" ...
## $ pheno_type : chr "type" "type" "type" "type" ...
analyses_tbl <- readRDS(file.path(projectpath, "analyses.rds"))
str(analyses_tbl)## 'data.frame': 26 obs. of 15 variables:
## $ pheno : chr "OF_distance_first4" "OF_distance" "OF_corner_pct" "OF_periphery_pct" ...
## $ longname : chr "OF_distance_first4" "OF_distance" "OF_corner_pct" "OF_periphery_pct" ...
## $ output : chr "OF_distance_first4" "OF_distance" "OF_corner_pct" "OF_periphery_pct" ...
## $ pheno_group: chr "group" "group" "group" "group" ...
## $ pheno_type : chr "type" "type" "type" "type" ...
## $ model : chr "normal" "normal" "normal" "normal" ...
## $ transf : chr "identity" "identity" "identity" "identity" ...
## $ offset : num 0 0 0 0 0 0 0 0 0 0 ...
## $ winsorize : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ sex : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ Cohort : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ Group : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ Subgroup : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ ngen : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ coat_color : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
The query routines embed the location and database information for
genes, variants, probabilities and mRNA data, which are generally very
large and only needed for a subset. All have first three arguments
(chr, start, end)
They are used in read_project and
read_query_rds, but we make their use explicit here for
encoding locations in query functions.
query_genes <- qtl2::create_gene_query_func(file.path(taxapath, "mouse_genes.sqlite"))
query_variants <- qtl2::create_variant_query_func(file.path(taxapath, "cc_variants.sqlite"))args("query_genes")## function (chr, start, end)
## NULL
str(query_genes(chr_id, scan_window[1], scan_window[2]))## 'data.frame': 11201 obs. of 15 variables:
## $ chr : chr "11" "11" "11" "11" ...
## $ source : chr "VEGA" "VEGA" "ENSEMBL" "ENSEMBL" ...
## $ type : chr "gene" "mRNA" "gene" "mRNA" ...
## $ start : num 96 96 96 96 96 ...
## $ stop : num 96 96 96 96 96 ...
## $ score : num NA NA NA NA NA NA NA NA NA NA ...
## $ strand : chr "-" "-" "-" "-" ...
## $ phase : num NA NA NA NA NA NA NA NA 2 NA ...
## $ ID : chr "VEGA:OTTMUSG00000001937" "VEGA:OTTMUST00000003877" "ENSEMBL:ENSMUSG00000013415" "ENSEMBL:ENSMUST00000013559" ...
## $ Name : chr NA NA NA NA ...
## $ Parent : chr NA "VEGA:OTTMUSG00000001937" NA "ENSEMBL:ENSMUSG00000013415" ...
## $ Dbxref : chr "MGI:MGI:1890357" "MGI:MGI:1890357,VEGA:OTTMUSG00000001937" "MGI:MGI:1890357" "MGI:MGI:1890357,ENSEMBL:ENSMUSG00000013415" ...
## $ mgiName: chr NA NA NA NA ...
## $ bioType: chr "KNOWN_protein_coding" NA "protein_coding" NA ...
## $ Alias : chr NA NA NA NA ...
args("query_variants")## function (chr, start, end)
## NULL
str(query_variants(chr_id, scan_window[1], scan_window[2]))## 'data.frame': 25588 obs. of 16 variables:
## $ snp_id : chr "rs218229099" "rs241313271" "rs27071476" "rs220108253" ...
## $ chr : chr "11" "11" "11" "11" ...
## $ pos : num 96 96 96 96 96 ...
## $ alleles : chr "A|T" "G|T" "A|G" "G|A" ...
## $ sdp : int 225 225 128 225 97 32 128 97 128 64 ...
## $ ensembl_gene: chr "ENSMUSG00000013415" "ENSMUSG00000013415" "ENSMUSG00000013415" "ENSMUSG00000013415" ...
## $ consequence : chr "intron_variant" "intron_variant" "intron_variant" "intron_variant" ...
## $ A_J : num 2 2 1 2 2 1 1 2 1 1 ...
## $ C57BL_6J : num 1 1 1 1 1 1 1 1 1 1 ...
## $ 129S1_SvImJ : num 1 1 1 1 1 1 1 1 1 1 ...
## $ NOD_ShiLtJ : num 1 1 1 1 1 1 1 1 1 1 ...
## $ NZO_HlLtJ : num 1 1 1 1 1 1 1 1 1 1 ...
## $ CAST_EiJ : num 2 2 1 2 2 2 1 2 1 1 ...
## $ PWK_PhJ : num 2 2 1 2 2 1 1 2 1 2 ...
## $ WSB_EiJ : num 2 2 2 2 1 1 2 1 2 1 ...
## $ type : chr "snp" "snp" "snp" "snp" ...
The query_probs assumes a probdir (default
genoprob), which contains the FST method files
probs_fstindex.rds (for genotype probabilities) and
aprobs_fstindex.rds (for allele probabilities). These are
FST master files, which point to individual FST databases by chromosome,
such as aprobs_1.fst and probs_1.fst.
query_probs <- qtl2pattern::create_probs_query_func(projectpath)args("query_probs")## function (chr = NULL, start = NULL, end = NULL, allele = TRUE,
## method = method_val, probdir = probdir_val)
## NULL
get("method_val", env = environment(query_probs))## [1] "fst"
get("probdir_val", env = environment(query_probs))## [1] "genoprob"
probs_obj <- query_probs(chr_id, scan_window[1], scan_window[2], allele = TRUE)
dim(probs_obj$probs[[1]])## [1] 261 8 9
range(probs_obj$map)## [1] 96.04762 97.95779
dim(query_probs(chr_id, scan_window[1], scan_window[2], allele = FALSE)$probs[[1]])## [1] 261 36 9
The query_mrna assumes a mrnadir (default
RNAseq), which should contain the following files:
annot.mrna.rds
annot.samples.rds
expr.mrna.fst
peaks.mrna.fst
query_mrna <- qtl2mediate::create_mrna_query_func(projectpath)The function has additional arguments local (use only
mRNA whose genes are in the region if TRUE) and
qtl (use only mRNA whose QTL peaks are in the region if
TRUE).
args("query_mrna")## function (chr = NULL, start = NULL, end = NULL, local = TRUE,
## qtl = FALSE, mrnadir = mrnadir_val)
## NULL
get("mrnadir_val", env = environment(query_mrna))## [1] "RNAseq"
There are no mRNA data for the Recla dataset, so this will not be used below. The following will create an error for these data.
str(query_mrna(chr_id, scan_window[1], scan_window[2]))
Phenotype mediation uses as mediators other phenotypes in the region
defined by chr_id and scan_window. This is all
bundled in the routine comediator_region, using
project_info and analyses_tbl settings. It
then calls pheno_region with peaks and
qtls.
comed_ls <- qtl2mediate::comediator_region(pheno_name, chr_id, scan_window,
covar, analyses_tbl, peaks,
qtls, pmap, pheno_data)med_ls <- qtl2mediate:::comediator_type(comed_ls, peaks, pheno_name, other)Need to either have query_mrna set up ahead of time (in
‘RDS’ format) or create here.
expr_ls <- get_expr_region(chr_id, scan_window[1], scan_window[2], covar, pmap,
drivers = qtls,
query_mrna = query_mrna)
Mediation test is performed in shinyMediate. Input includes
med_ls (from other phenotypes) or expr_ls
(from mRNA expression data). Genetic cross data has kinship information,
which is important to adjust for genetic correlation among
individuals.
kinship <- readRDS(file.path(projectpath, "kinship.rds"))[[chr_id]]The position for target mediation is by default the center of the map
in this region. However, it might be improved by looking at peak for
pheno_name.
pos_Mbp <- mean(range(probs_obj$map))peak_mar <- qtl2::find_marker(probs_obj$map, chr_id, pos_Mbp)
geno_max <- subset(probs_obj$probs, chr = chr_id, mar = peak_mar)[[1]][,,1]cov_tar <- qtl2mediate:::get_covar(
covar,
dplyr::filter(analyses_tbl, pheno == pheno_name))
driver_med <- probs_obj$probs[[chr_id]]phe_mx <- pheno_data[, pheno_name, drop = FALSE]med_scan <- intermediate::mediation_scan(
target = phe_mx,
mediator = med_ls[[1]],
driver = geno_max,
annotation = med_ls[[2]],
covar = cov_tar,
kinship = kinship,
fitFunction = qtl2mediate::fitQtl2)summary(med_scan, minimal = TRUE)## chr pos LR
## LD_distance_light 11 97.02702 3.188115
## LD_transitions 11 96.68386 8.798000
## VC_bottom_distance 11 97.02702 9.884018
## VC_top_distance_first4 11 96.17147 10.321125
## OF_distance 11 97.02702 11.709825
## VC_bottom_transitions 11 97.29928 12.331006
## OF_distance_first4 11 97.02702 12.596732
ggplot2::autoplot(med_scan, size = 3)med_test <- qtl2mediate::mediation_test_qtl2(
target = phe_mx,
mediator = med_ls[[1]],
annotation = med_ls[[2]],
covar_tar = cov_tar,
covar_med = med_ls$covar,
genoprobs = probs_obj$probs,
map = probs_obj$map,
chr = chr_id,
pos = pos_Mbp,
kinship = kinship)summary(med_test)## # A tibble: 7 × 30
## id chr pos triad pvalue alt_t…¹ undec…² LR_me…³ LR_me…⁴ LR_CMST df
## <chr> <fct> <dbl> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LD_dis… 11 97.0 caus… 0 reacti… 0 3.19 15.1 141. 8
## 2 LD_tra… 11 96.7 caus… 0 reacti… 0 8.8 13 31.5 8
## 3 OF_dis… 11 97.0 caus… 0 reacti… 0 11.7 10.6 21.8 8
## 4 OF_dis… 11 97.0 caus… 0 reacti… 0 12.6 11.5 20.0 8
## 5 VC_bot… 11 97.0 caus… 0 reacti… 0 11.1 11.7 28.6 8
## 6 VC_bot… 11 97.3 caus… 0 reacti… 0 13.5 8.5 20.8 8
## 7 VC_top… 11 96.2 caus… 0 reacti… 0 11.5 11.1 18.8 8
## # … with 19 more variables: IC <dbl>, lod <dbl>, longname <chr>, output <chr>,
## # pheno_group <chr>, biotype <chr>, model <chr>, transf <chr>, offset <dbl>,
## # winsorize <lgl>, sex <lgl>, Cohort <lgl>, Group <lgl>, Subgroup <lgl>,
## # ngen <lgl>, coat_color <lgl>, qtl_ct <int>, info <chr>, local <lgl>, and
## # abbreviated variable names ¹alt_test, ²undecided, ³LR_mediation,
## # ⁴LR_mediator
ggplot2::autoplot(med_test, size = 3, alpha = 1)Based on shinyTriad in qtl2shiny. Basic idea is to find traits with strongest additive allele and/or SNP signal in the region based on scans. Then uses these to identify the best SDP(s) for further investigation with triads.
snpinfo <- query_variants(chr_id, scan_window[1], scan_window[2])
snpprobs <- qtl2mediate:::get_snpprobs(chr_id, peak_Mbp, window_Mbp,
pheno_name,
probs_obj$probs,
probs_obj$map,
snpinfo)snp_scan_obj <- qtl2mediate:::scan1covar(pheno_data, cov_tar,
snpprobs$snpprobs, kinship,
analyses_tbl,
sex_type = "all")qtl2ggplot::ggplot_scan1(snp_scan_obj, snpprobs$snpinfo,
lodcolumn = pheno_name,
patterns = "hilit", drop_hilit = 3,
haplos = snpprobs$haplos,
main = pheno_name)Want to use qtl2::top_snps instead of qtl2pattern::top_snps_pattern,
but then have to do the work. Put this in triad_helper() in
’medation_triad_qtl2.R`.
(tmp_snps_1 <-
dplyr::arrange(
qtl2::top_snps(snp_scan_obj,
snpprobs$snpinfo,
lodcolumn = pheno_name,
show_all_snps = FALSE),
dplyr::desc(.data$lod)))## snp_id chr pos alleles sdp ensembl_gene consequence A_J
## 1 rs257161856 11 97.00586 G|A 116 intergenic_variant 1
## 2 rs27023987 11 97.00335 C|G/T 139 intergenic_variant 3
## C57BL_6J 129S1_SvImJ NOD_ShiLtJ NZO_HlLtJ CAST_EiJ PWK_PhJ WSB_EiJ type index
## 1 1 2 1 2 2 2 1 snp 12491
## 2 2 1 2 1 1 1 2 snp 12464
## interval on_map lod
## 1 3 FALSE 4.884181
## 2 3 FALSE 4.884181
qtl2mediate:::sdp_to_pattern(tmp_snps_1$sdp[1:2], haplos = snpprobs$haplos)## [1] "ABDH:CEFG" "ABDH:CEFG"
Not quite what we want here. Want top SDP and top peaks; 1 or more.
snp_id_sdp <- dplyr::filter(snpprobs$snpinfo, sdp %in% tmp_snps_1$sdp[1:2])$snp_id
m <- match(snp_id_sdp, rownames(snp_scan_obj), nomatch = 0)
tmp <- as.data.frame(t(snp_scan_obj[m,]))
dplyr::arrange(tmp, dplyr::desc(.data[[names(tmp)[1]]]))## rs27023987 rs257161856
## LD_light_pct 4.884180816 4.884180816
## LD_distance_light 4.290708137 4.290708137
## LD_transitions 3.536163536 3.536163536
## VC_top_distance_first4 3.533595886 3.533595886
## VC_top_distance 2.592797951 2.592797951
## OF_immobile_pct 1.932726913 1.932726913
## VC_top_velocity 1.486847315 1.486847315
## OF_distance 1.042637899 1.042637899
## VC_bottom_distance 1.006557071 1.006557071
## TS_frequency_climbing 0.952787113 0.952787113
## HP_latency 0.717622729 0.717622729
## OF_corner_pct 0.702872663 0.702872663
## VC_bottom_transitions 0.584050714 0.584050714
## OF_distance_first4 0.550941369 0.550941369
## VC_bottom_velocity 0.503672483 0.503672483
## TS_time_immobile 0.441407270 0.441407270
## LD_light_first4 0.406734913 0.406734913
## OF_center_pct 0.392502502 0.392502502
## bw 0.359302531 0.359302531
## VC_bottom_distance_first4 0.292490619 0.292490619
## OF_periphery_pct 0.224770141 0.224770141
## VC_top_time_pct 0.147379093 0.147379093
## VC_bottom_time_pct 0.143988264 0.143988264
## TS_latency_immobile 0.143689933 0.143689933
## VC_bottom_time_first4 0.008482029 0.008482029
## VC_top_time_first4 0.003921065 0.003921065
qtl2ggplot::ggplot_scan1(snp_scan_obj, snpprobs$snpinfo,
lodcolumn = pheno_name,
patterns = "hilit", drop_hilit = 3,
haplos = snpprobs$haplos)qtl2ggplot::ggplot_scan1(snp_scan_obj, snpprobs$snpinfo,
lodcolumn = match(c("LD_light_pct",
"LD_transitions",
"LD_distance_light"),
colnames(snp_scan_obj)),
#seq_len(ncol(snp_scan_obj))[1:2],
patterns = "hilit", drop_hilit = 3,
facet = "pheno",
haplos = snpprobs$haplos)top_snps_tbl <- qtl2pattern::top_snps_pattern(
snp_scan_obj,
snpprobs$snpinfo)summary(top_snps_tbl)## # A tibble: 36 × 8
## pheno min_pos max_pos max_lod min_lod sdp pattern snp_id
## <chr> <dbl> <dbl> <dbl> <dbl> <int> <chr> <chr>
## 1 LD_light_pct 97.0 97.0 4.88 4.88 116 ABDH:CEFG rs257…
## 2 LD_light_pct 97.0 97.0 4.88 4.88 139 ABDH:CEFG rs270…
## 3 LD_transitions 96.7 96.8 4.46 3.42 96 ABCDEH:FG 1819 …
## 4 VC_top_distance_first4 96.7 96.7 4.36 3.82 138 BDH:ACEFG 19 SN…
## 5 LD_distance_light 97.0 97.0 4.29 4.29 116 ABDH:CEFG rs257…
## 6 LD_distance_light 97.0 97.0 4.29 4.29 139 ABDH:CEFG rs270…
## 7 LD_distance_light 97.0 97.0 3.86 3.40 101 BDEH:ACFG 8 SNPs
## 8 VC_top_distance_first4 96.7 96.7 3.82 3.82 117 BDH:ACEFG rs249…
## 9 VC_top_distance_first4 96.6 96.7 3.74 3.74 101 BDEH:ACFG 6 SNPs
## 10 LD_transitions 97.0 97.0 3.54 3.54 116 ABDH:CEFG rs257…
## # … with 26 more rows
patterns <- summary(top_snps_tbl)
haplos <- LETTERS[1:8]
triad <- unique(med_test$best$triad)
# Pick top trait with causal call.
med_name <- dplyr::filter(med_test$best, .data$triad == "causal")[["id"]][1]
pattern <- dplyr::filter(patterns, .data$pheno == med_name)$pattern[1]
sdps <- unique(dplyr::filter(patterns, .data$pheno == med_name)$sdp)
sdp <- sdps[qtl2mediate:::sdp_to_pattern(sdps, haplos) == pattern]
id <- med_ls[[2]]$id[med_ls[[2]][["id"]] == med_name]med_triad <- intermediate::mediation_triad(target = phe_mx,
mediator = med_ls[[1]][, id, drop = FALSE],
driver = geno_max,
covar_tar = cov_tar,
covar_med = med_ls$covar,
kinship = kinship,
sdp = sdp, allele = TRUE,
fitFunction = qtl2mediate::fitQtl2)Look carefully at code in last section, and at
triad_sdp. Is the goal to find the sdp for the best triad
across mediators or is it to just look at preselected mediator?
sdp <- qtl2mediate::triad_sdp(chr_id, pos_Mbp,
scan_window,
pheno_name, pheno_data = pheno_data,
covar_tar = cov_tar,
probs_obj$probs, probs_obj$map,
analyses_tbl,
kinship = kinship)Need to customize to particular mediator chosen.
med_triad <- qtl2mediate::mediation_triad_qtl2(
target = phe_mx,
mediator = med_ls[[1]][, id, drop = FALSE],
annotation = med_ls[[2]],
covar_tar = cov_tar,
covar_med = med_ls$covar,
genoprobs = probs_obj$probs,
map = probs_obj$map,
chr = chr_id,
pos = pos_Mbp,
kinship = kinship,
sdp)summary(med_triad)## Model fit with driver columns
##
## |term | estimate| std.error| statistic| p.value|
## |:----------------|-----------:|---------:|----------:|---------:|
## |mediator | 0.0062485| 0.0026272| 2.3784080| 0.0181488|
## |SexFemale | 35.5354522| 5.6471836| 6.2925973| 0.0000000|
## |SexMale | 33.7103908| 6.1341315| 5.4955442| 0.0000001|
## |A | 2.0962236| 6.1885618| 0.3387255| 0.7351043|
## |B | 0.7602089| 6.3443328| 0.1198249| 0.9047193|
## |C | -15.6357536| 6.1863361| -2.5274659| 0.0121137|
## |D | -1.4510211| 5.9428114| -0.2441641| 0.8073064|
## |E | -4.7513767| 5.9949686| -0.7925607| 0.4287944|
## |F | -14.8780078| 5.4475045| -2.7311603| 0.0067662|
## |G | -14.7548869| 5.6738148| -2.6005232| 0.0098693|
## |H | 0.0000000| 0.0000000| 0.0000000| 0.0000000|
## |mediator:SexMale | 0.0028075| 0.0038471| 0.7297690| 0.4662229|
##
## Model fit with driver groups
##
## |term | estimate| std.error| statistic| p.value|
## |:----------------|-----------:|----------:|----------:|---------:|
## |mediator | 0.0054653| 0.0028453| 1.9208470| 0.0560495|
## |SexFemale | 42.7192755| 8.5793636| 4.9793059| 0.0000013|
## |SexMale | 40.3536882| 8.6381221| 4.6715811| 0.0000052|
## |groupAB | 0.8476695| 10.6113354| 0.0798834| 0.9364029|
## |groupAC | -19.1925195| 8.8554365| -2.1673149| 0.0312895|
## |groupAD | -21.1649362| 10.6412234| -1.9889570| 0.0479506|
## |groupAE | -6.5846890| 16.9331714| -0.3888633| 0.6977551|
## |groupAF | -19.0213010| 8.8830746| -2.1412970| 0.0333546|
## |groupAG | -4.6151459| 8.5793530| -0.5379364| 0.5911672|
## |groupAH | -5.1486070| 8.7570079| -0.5879413| 0.5571774|
## |groupBB | 2.9607322| 16.7707948| 0.1765410| 0.8600321|
## |groupBC | -5.1323361| 9.3957988| -0.5462373| 0.5854588|
## |groupBD | -3.7704978| 16.9150829| -0.2229074| 0.8238153|
## |groupBE | -16.7587658| 10.5788947| -1.5841698| 0.1145980|
## |groupBF | -14.7930949| 8.2749856| -1.7876883| 0.0752090|
## |groupBG | -16.7602294| 9.0182555| -1.8584780| 0.0644429|
## |groupBH | -5.2897217| 9.7078816| -0.5448894| 0.5863840|
## |groupCC | 7.2480651| 16.7691641| 0.4322258| 0.6660026|
## |groupCD | -11.8523945| 10.0588926| -1.1783001| 0.2399557|
## |groupCE | -15.0919469| 8.8576031| -1.7038409| 0.0898292|
## |groupCF | -29.5018343| 9.2084014| -3.2037954| 0.0015584|
## |groupCG | -28.0035295| 8.9753249| -3.1200575| 0.0020514|
## |groupCH | -7.3006264| 9.4034348| -0.7763787| 0.4383631|
## |groupDD | -7.4664462| 9.6836985| -0.7710325| 0.4415187|
## |groupDE | -2.9445950| 8.6502080| -0.3404074| 0.7338760|
## |groupDF | -17.6802963| 8.5688101| -2.0633315| 0.0402602|
## |groupDG | -22.6306973| 10.0363535| -2.2548725| 0.0251303|
## |groupDH | -2.1399215| 9.3840554| -0.2280380| 0.8198294|
## |groupEE | -11.2954916| 11.4532466| -0.9862262| 0.3251105|
## |groupEF | -19.4211281| 9.2044763| -2.1099656| 0.0359964|
## |groupEG | -19.2922730| 9.2171838| -2.0930767| 0.0374936|
## |groupEH | -11.5139640| 16.7588549| -0.6870376| 0.4927859|
## |groupFF | -12.8569264| 9.1678473| -1.4023932| 0.1622138|
## |groupFG | -15.7006508| 8.2758211| -1.8971714| 0.0591210|
## |groupFH | -15.0044632| 8.4325119| -1.7793587| 0.0765678|
## |groupGG | -29.5473583| 11.4551739| -2.5793898| 0.0105518|
## |groupGH | -4.4104300| 10.6128505| -0.4155745| 0.6781283|
## |groupHH | -23.3427702| 11.6163964| -2.0094674| 0.0457148|
## |mediator:SexMale | 0.0035596| 0.0040655| 0.8755636| 0.3822266|
intermediate::ggplot_mediation_triad(med_triad)intermediate::ggplot_mediation_triad(med_triad, fitlines = "parallel")plotly::ggplotly(intermediate::ggplot_mediation_triad(med_triad, fitlines = "parallel"))